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The QCD phase diagram at nonzero quark density 
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Abstract: We determine the phase diagram of QCD on the \x — T plane for small to 
moderate chemical potentials. Two transition lines are defined with two quantities, the 
chiral condensate and the strange quark number susceptibility. The calculations are carried 
out on Nt = 6, 8 and 10 lattices generated with a Symanzik improved gauge and stout-link 
improved 2+1 flavor staggered fermion action using physical quark masses. After carrying 
out the continuum extrapolation we find that both quantities result in a similar curvature 
of the transition line. Furthermore, our results indicate that in leading order the width of 
the transition region remains essentially the same as the chemical potential is increased. 
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1. Introduction 

The understanding of the phase diagram of QCD is of utmost importance and has attracted 
much attention, both experimental and theoretical. Experimental results are coming from 
cosmology and heavy ion collisions. Recently, in a collision of gold nuclei at the Relativis- 
tic Heavy Ion Collider (RHIC), a temperature beyond 200 MeV was reached [1], which 
indicates that the quark-gluon plasma has been created. Furthermore, the density of the 
system can be varied by tuning the center of mass energy ^snn- While most of the on- 
going experiments like those at LHC or RHIC concentrate on achieving very high energies 
and thus small chemical potentials, there are projects that aim for regions of the phase dia- 
gram with larger densities (RHIC II, Facility for Antiproton and Ion Research (FAIR)). In 
these latter experiments an important objective is to identify the critical endpoint, e.g., by 
searching for critical opalescence. Designing these next generation experiments can benefit 
greatly from developing theoretical understanding of the phase diagram. 

Our theoretical knowledge about the phase diagram of QCD is mostly limited to the 
zero chemical potential (p = 0) axis and obtained by the use of lattice QCD. The main 
reason that full results for fi > are not available is the infamous sign problem which spoils 
any lattice technique based on importance sampling. There are various scenarios for the 
fi > region of the phase diagram, among which two are illustrated in Figure 1. 

The transition at fi = is a crossover [2] and we expect that the transition temperature 
decreases as we increase fx. Besides the actual value of the curvature of the transition line a 
particularly interesting question is whether the transition becomes weaker or stronger as fi 
grows. A strengthening of the transition could lead to the existence of a critical point, where 
the crossover transforms into a true phase transition (see left side of Figure 1). Another 
possibility is that the transition weakens with increasing \x (see right side of Figure 1). The 
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Figure 1: Two possible scenarios for the QCD phase diagram on the \x — T plane, defined using 
a given observable. The left panel shows a phase diagram with a transition growing stronger and 
possibly even turning into a real, first-order phase transition at a critical endpoint. The right panel 
on the other hand corresponds to a scenario with a weakening transition and no critical endpoint. 
The paths corresponding to systems describing the early Universe and a heavy ion collision are also 
shown by the arrows. Note that different observables may lead to different scenarios. 



existence of the critical point would not be ruled out by such a scenario but would require 
non-monotonic behavior [3]. 

At zero chemical potential lattice calculations provide reliable and accurate results [4— 
9]. Much more difficult is the situation at nonzero chemical potential. Simulations at 
non-vanishing chemical potential are burdened by the sign (complex action) problem: the 
fermion determinant here becomes complex, and as a result makes Monte-Carlo methods 
based on importance sampling impossible. Recently several methods were developed to 
access the region of small chemical potentials. They are all based on simulations at zero or 
purely imaginary chemical potentials where the sign problem is absent. The first possibility 
is the reweighting of the generated configurations [10-14]. The weight factors can also be 
approximated by a Taylor expansion in \i [15-19]. Further possibilities are an analytic 
continuation from imaginary \i [20-27], or using the canonical ensemble [28-30]. The 
above studies were carried out on coarse lattices and in most cases with non-physical quark 
masses. We emphasize that to have a full result, the use of physical quark masses and a 
reliable continuum extrapolation are essential. In this paper we determine the transition 
temperature T c (/x) as a function of the chemical potential through a Taylor-expansion 
technique. The first term of this expansion is zero due to the symmetry Z{\i) = Z(—fi) 
of the partition function. Therefore the first nonvanishing contribution comes from the 
second order, which is related to the curvature of the transition line. 



2. Definition of the curvature 

Let us parameterize the transition line in the vicinity of the vertical fx = axis as 

T c ( fJ , 2 ) = T c (l-K-(i 2 /T?) (2.1) 
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with T c being short for T c (0). This implies that the curvature can be written as 



where (and also in the following) fi refers to the baryonic chemical potential (fi = fiB = 
3/J"u,d), where fj, Uj d is the quark chemical potential assigned to the light quarks. Thus one 
has to measure T c as a function of \x for small chemical potentials. To this end we use a 
definition of T c which is most suitable for determining the curvature. 

Let us consider a quantity <j)(T,y?) that is monotonic in T in the transition region, 
and fulfills the following constraints: 



lim -— ^(T, z" 2 ) = 0, lim ^(f>(T, fi 2 ) = (2.3) 

T^O <9/i 2 T^oo d/J, 2 



that is to say, <f> does not depend on the chemical potential in the limiting cases T — > and 
T — > co. For any fixed /x we can define a transition temperature T c {fi 2 ) as the temperature 
at which 4>(T,/j, 2 ) takes the predefined constant value C: 

HT,^)\ T=W) = C. (2.4) 

We will choose a C that corresponds to the inflection point of <f>(T,0). (Note that T c can 
also be defined as the location of the maximum or inflection point of some observable. At 
non-zero \x this turns out to be somewhat less advantageous since a fitting of the reweighted 
data is required.) 

Now let us determine the curvature using this definition of T c (fj, 2 ). The total derivative 
of the observable (f>(T, ju 2 ) may be written as 

d^ = (d4>/dT)\^ =0 ■ at + (d<f>/d(v 2 )) | M=0 • d^ 2 (2.5) 

Along the T c {n 2 ) line, 4> is constant by definition, thus d(j) = 0. One obtains 

dr c f d<f>\ /(d<f>\ 

T=T C 
H=0 



d^ 2 \d/J, 



t=t c I \ dT 



R(T) 

Thus, for every C we can define a curvature. Since the T C {C) function is invertible for the 
whole C range, we can also write (2.6) as a function of temperature, R(T). 

The function R(T) is related to the distance that the 4>(T) curve shifts along the T 
axis as the chemical potential is varied. Given <p(T) and R(T) at zero chemical potential, 
the shift for non-zero fi at leading order is R(T) ■ fi 2 (the curve moves to the left if R(T) is 
negative and to the right otherwise). This behavior is illustrated in Figure 2. Using R(T) 
we can define a temperature dependent curvature according to (2.2) as k(T) = —T c • R(T). 
The meaning of k(T) is again simple: it gives the curvature of the 4> = const, curve which 
starts from T at fi = 0. 

We use the value of k(T) at T = T c to define the curvature for a given observable. The 
shape of the k(T) function also has important consequences. The slope of n(T) around 
T c is related to the width of the transition as follows: if the slope is zero, i.e. n(T) is 
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Figure 2: Illustration of the behavior of the observable at /Lt = and /z > 0. The quantity T c (/i 2 ) 
is defined as the temperature where 4>(T,n 2 ) crosses a constant value C. For fj, > 0, each point of 
the </>(T, 0) curve shifts in T by R(T) ■ /i 2 (see definition in text). 



constant around T c , then all points shift the same amount along the T axis when a small 
chemical potential is switched on. This means that to leading order in \i the shape of 
the 4>(T) function (and thus, the width of the transition) does not change. If the slope is 
positive, then points with larger T shift more than the ones with smaller T resulting in a 
compression of points, i.e. a narrower transition. Similarly, a negative slope indicates a 
broadening of the transition. 

All in all, the expression dn/dT is therefore related to the relative change in the width 
W([a) of the transition as the chemical potential increases: 



1 dW 1 8k 



W d(/i 2 ) T c dT 



T=T C 



where we assume that W is proportional to the inverse slope of the quantity in question: 

w~\{dcj>idT\ T=Ta )- x \. 

The two observables we use are the renormalized chiral condensate 4> = (ipip r ) an d the 
normalized strange quark number susceptibility <p = {Xs/T 2 }. As we show in section 4, 
both satisfy the constraints listed in (2.3). The derivative d(/)/dT is determined numerically, 
using the fi = data as a function of the temperature. In order to calculate the derivative 
d(p/d(fi 2 ) we need to measure more complicated operators; the technique for computing 
these is detailed in the next section. 



3. Determination of the Taylor-coefficients 

Let us consider the partition function of the staggered lattice formulation for Nf fermion 
flavors in its usual form 

Z = J VUe- s ° iu \detM) N f^ (3.1) 

and denote the derivative with respect to [i u ^ by '. The derivatives of Z are easily calculated 
to be (log Z)' = (n Ui d) and (log Z)" = (r? u d + n' u d ) — (n Ui d) 2 , where the light quark number 
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density n u ^ and its derivative with respect to jjL u ^ are given by the following combinations: 

n u4 = -fTr {M- l M') 
Nf 

n ' u d = _2_Tr {M~ l M" - M^M'M^M') 

Using these definitions the second derivative of any (possibly explicitly //-dependent) ob- 
servable can be straightforwardly determined. For the renormalized chiral condensate and 
the strange quark number susceptibility (see definition in section 4) one obtains: 

= Wr«d + <d)> " Wr)« d + <d> + {^' r nu,d + Wt) (3-2) 

= (X«("u,d + n 'u,d)) ~ (Xs){nl td + n' u4 ) - 2{n s n u4 ) 2 (3.3) 

where n s is the strange quark number density, defined similarly as n u ^. Note that the ad- 
ditive renormalization of V>Vv ( see section 4.2) does not influence the derivative in question. 

For the chiral condensate - being a //-dependent operator - the derivatives ipip' r and 
V'Vr of this explicit dependence are also present in (3.2). These terms were calculated 
numerically, using a purely imaginary chemical potential A/tj. The value of A//j was varied 
in the range 0.01 . . . 0.0005, and it was checked that the finite differences converge fast 
enough to the A//j — > values and the error coming from this approximation is negligible 
compared to statistical errors. Taking into account these considerations A//j = 0.001 was 
used. 



d 2 (Xs) 



4. The //-dependence of the observables 

We calculated the curvature of the transition line using the strange quark number suscep- 
tibility and the chiral condensate. The details of their renormalization and behavior are 
explained in this section. 

4.1 The strange quark number susceptibility 

The strange quark number susceptibility is defined as 

T dHogZ 

This observable needs no renormalization, since it is connected to a conserved current. It 
is useful to study the combination (x s /T 2 ), since it obeys the conditions of (2.3). It is easy 
to see that at T = one gets (xs/T 2 ) = and at T — > oo the normalized quark number 
susceptibility (x s /T 2 ) reaches its fi Uj d independent Stefan-Boltzmann limit of 1. 

4.2 The chiral condensate 

The chiral condensate can also be expressed as a derivative of the partition function: 

<*> - ™ (-) 
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The renormalization of tpt/j is a more delicate issue as compared to the situation with \s- 
The free energy (log Z) contains additive divergences in the cutoff. In order to carry out the 
proper renormalization of the condensate, these additive divergences have to be eliminated 
- this is done by subtracting the T = contribution. 

The multiplicative divergence due to the derivative with respect to the mass can be 
eliminated with a multiplication by the bare quark mass. Then, in order to have a dimen- 
sionless combination 1 , the whole expression can be divided by the fourth power of some 
dimensionful mass scale, Q 4 : 



0))-m-ij 



(4.3) 



This way no divergent contributions remain: this is a meaningful quantity to study in the 
continuum limit. 2 In this work we use the T = pion mass for the Q normalization scale. 

The final condition that has to be satisfied is that {ipipr) should be independent of 
fi at T = and T — >• oo. At T = the partition function is independent of ji as long 
as n is smaller than a /j, c critical value (the approximate baryon mass) and no baryons 
can be created from the vacuum. Only for [i > [i c does the partition function have a 
non-trivial \x dependence. Therefore all derivatives of Z (thus (tptpj.)) are independent of 
fj, for fj, < /j, c . The chemical potential regime covered in this paper lies in this region. In 
the Stefan-Boltzmann limit (T — > oo) the (i independence is only satisfied in the sense 



(tptpr) 1 d/dfj 2 (ip^ r 



0. Figure 4 demonstrates, however, that for temperatures above 



the transition region ipij) r is already practically independent of [i. 

Figure 3 illustrates the behavior of (xs/T 2 ) and (V>Vv) a s a function of the temperature, 
determined on N t = 10 lattices. 
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Figure 3: The strange susceptibility and the renormalized chiral condensate as functions of the 
temperature at jj. — 0. 



1 Note that a division by T 4 which would also render the condensate dimensionless, changes the temper- 
ature dependence and would lead to a non-monotonic T dependence, which would be disadvantageous in 
the present context. 

2 Note that this renormalization procedure leads to a somewhat unusual chiral condensate which vanishes 
at T — and reaches a negative value at T — > oo. A more conventional condensate which is positive at 
T = and goes to zero at large temperatures can be obtained by a constant shift which is irrelevant for 
our present study. 
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5. Simulation setup 



We used a Symanzik improved gauge and stout-link improved staggered fermionic lattice 
action in order to reduce taste violation [5]. The configurations were generated with an 
exact RHMC algorithm [31]. We determined the line of constant physics (LCP) using 
physical masses for the light quarks m u ^ as well as for the strange quark m s . The LCP 
was fixed by setting the ratios mx/fx and mx/inn to their physical values. We used three 
different lattice spacings Nt = 6, 8, 10 and aspect ratios N s /Nt of 4 and 3. The scale was 
fixed by fx and its unambiguity was checked by calculating uik*, fn and r$. The random 
noise estimator method was used to measure the operators detailed in section 3. We used 
80 random vectors so that the error coming from the method and the statistical error are 
of the same extent. The details of the simulation setup can be found elsewhere [5,32]. 
We used the gauge ensembles generated for a \i = study [5]. We also generated extra 
configurations for N t = 8 and 10. The number of trajectories at each fj, N s and Nt is 
summarized in table 1. The autocorrelation times were below 10 trajectories in all cases. 
After confirming the absence of thermalization effects, we measured observables on every 
fifth trajectory. The measurements were performed on clusters equipped with graphics 
cards [33]. 



Nf x N t 
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#of trajecs. 


Nf x N t 
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#of trajecs. 
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#of trajecs. 




3.4500 


1750 




3.6000 


1800 




3.6500 


800 




3.4950 


2500 




3.6250 


2100 




3.6750 


3350 




3.5100 


5200 
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3.7000 


800 




3.5250 


5350 




3.6500 


5000 




3.7125 


850 


24 3 x 6 


3.5400 


5450 




3.6625 


3150 




3.7250 


800 


3.5550 


3400 


24 3 x 8 


3.6750 
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3.5700 
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2100 


3.7750 
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3.6450 
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3.7160 


3050 




3.7875 


4200 


18 3 x 6 


3.5550 


4550 




3.7400 


2850 




3.8000 


1600 














3.8125 


2000 














3.8250 


4850 














3.8375 


4150 














3.8550 


5950 



Table 1: Number of trajectories for various lattice geometries. 



6. Results 

First we checked finite size effects by comparing our results at /3 = 3.555 obtained on 24 3 x 6 
and on 18 3 x 6 lattices. This value of f3 corresponds to about 155 MeV, i.e. is near the 
pseudocritical temperature. The larger box is of physical size ~ 5 fm. We observe a good 
agreement as the results for d4>/d{^ 2 ) agree within statistical errors for both the chiral 
condensate <f> = (tpip r ) and the strange quark number susceptibility <p = (x s /T 2 ). Figure 4 
shows our Nt = 6 results for N s = 24 and N s = 18. Thus we conclude that finite size errors 
can be neglected at the present statistical accuracy. 

Since the actual shape of the k(T) function is unknown we carry out a Taylor expansion 
around T c in the t = (T — T c )/T c dimensionless variable: 

k(T) = k(T c ) + c • t + ci • t 2 (6.1) 
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Figure 4: The derivative of our observables with respect to ^ measured on N t = 6 lattices. The 
N s = 24 results (blue points) are checked at one temperature by N s = 18 (red point). In the case 
of both observables a good agreement is observed, which indicates that finite size effects are small 
as compared to statistical errors. 



For each lattice spacing (i.e. each Nt) we have several simulation points, corresponding to 
different temperatures. In order to fit all of our points at once, we allow a lattice spacing 
dependence for the constant and linear terms (having a lattice spacing dependence of the 
quadratic term is also possible, but it does not improve the quality of the fits). Therefore 
we fit all of our simulation points with the following function 

k(T; N t ) = k(T c ; cont) + c • t + ci • t 2 + c 2 /iV f 2 + c 3 • t/N? (6.2) 

with fit parameters k(T c ; cont), cq, c\, Ci and C3. The independent data points as well as 
the fitted curves (for each JVj and in the continuum) are shown in Figure 5. 

The x 2 /d.o.f. values of the two fits are 1.19 and 1.29, respectively, indicating good fit 
qualities. The continuum curvatures are given by the k(T c ; cont) fit parameter, while the 
relative change in the width of the transition can be read off from — cq. Our final results 
are 

k (xs/t*) = 0.0089(14), = 0.0066(20) 

AW/W^ Xs/t2) = 0.033(16), AW/W^ = 0.030(18) 

The results obtained from the two quantities are consistent with each other. Using the n 
values we can give the transition lines defined by any of the observables as 

T e (p) = r qM=0 [1 - * ■ ^/Tc-^o] (6-3) 

The results for AW/W also suggest that the transition remains a weak crossover with 
essentially constant strength for small to moderate chemical potentials. Actually, there 
is a slight increase in the width of the transition determined from both quantities. This 
effect is, however, very weak: the width only changes by a few percent up to [i ~ T c . This 
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Figure 5: The curvature n(T) (see definition in text) determined using the strange quark number 
susceptibility (left) and the renormalized chiral condensate (right), respectively. A result of the 
combined fit (described in the text) is shown by the gray band. The fit results for the individual 
N t — 6, 8, 10 lattices are shown by the red, blue and green curves. The width of the gray band 
corresponds to the statistical uncertainty of the fit. 



finding is consistent with previous results in the literature. In Ref. [13] the imaginary parts 
of the Lee- Yang zeroes of the partition function were studied with exact reweighting. At 
sufficiently high fi the Lee- Yang zeroes approached the real axis, thus leading to a critical 
point. However, around \x = an opposite effect can be observed, the Lee- Yang zeroes 
slightly move away from the real axis, indicating a weakening of the transition. In Refs [25] 
(3 flavors) and [27] (2+1 flavors) the critical surface in the quark masses - chemical potential 
space was studied and the curvature of the surface suggests a weakening of the transition as 
the chemical potential is increased. Since all these leading order results predict a weakening 
of the transition for real chemical potentials (i.e. fi 2 > 0), a strengthening is expected in 
the imaginary direction (fi 2 < 0). It is interesting to note that these leading order results 
are of the same order of magnitude in the sense that using an extreme extrapolation they 
all lead to a critical point for an imaginary chemical potential fj, q = i ■ m in the range 
Hj/T « 1 - 3. 

The validity range of our result is difficult to estimate from the present study alone. A 
conservative estimate for the limit where the result obtained through the Taylor-expansion 
is still reliable is [Xu,d ~ T c i.e. where the expansion parameter exceeds unity. In the 
baryonic chemical potential this corresponds to about 500 MeV. Beyond this limit higher- 
order corrections are by all means expected to be important. Earlier experience with the 
exact reweighting method [13] also shows that the leading-order quadratic behavior of 
the T c {jj) function dominates upto the above mentioned limit in the baryonic chemical 
potential. To investigate whether higher order terms may lead to a critical point one must 
carry out a similar analysis with full reweighting, beyond the reach of present computational 
resources. 

Our final result is shown in Figure 6. The crossover region's extent changes little as the 
chemical potential increases, and within it two definitions give different curves for T c (fj,). 
It is useful to compare the whole picture to the freeze-out curve [34] which summarizes 
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experimental results on the {T, //} points where hadronization of the quark-gluon plasma 
was observed. This curve is expected to lie in the interior of the crossover region, as is 
indicated by our results as well. 
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Figure 6: The crossover transition between the 'cold' and 'hot' phases is represented by the 
coloured area (blue and red correspond to the transition regions obtained from the chiral condensate 
and the strange susceptiblity, respectively). The lower solid band shows the result for T c (/x) defined 
through the chiral condensate and the upper one through the strange susceptibility. The width 
of the bands represent the statistical uncertainty of T c (jx) for the given \i coming from the error 
of the curvature k for both observables. The dashed line is the freeze-out curve from heavy ion 
experiments [34]. Also indicated are with different symbols the individual measurements of the 
chemical freeze-out from RHIC, SPS (Super Proton Synchrotron) and AGS (Alternating Gradient 
Synchrotron), respectively. The center of mass energies ^/snn for each are shown in the legend. 
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